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ABSTRACT Bacteria living on the aerial parts of plants (the phyllosphere) are globally abundant and ecologically significant com- 
munities and can have significant effects on their plant hosts. Despite their importance, little is known about the ecological pro- 
cesses that drive phyllosphere dynamics. Here, we describe the development of phyllosphere bacterial communities over time on 
the model plant Arabidopsis thaliana in a controlled greenhouse environment. We used a large number of repUcate plants to 
identify repeatable dynamics in phyllosphere community assembly and reconstructed assembly history by measuring the com- 
position of the airborne community immigrating to plant leaves. We used more than 260,000 sequences from the v5v6 hyper- 
variable region of the 16S rRNA gene to characterize bacterial community structure on 32 plant and 21 air samples over 73 days. 
We observed strong, reproducible successional dynamics: phyllosphere communities initially mirrored airborne communities 
and subsequently converged to a distinct community composition. While the presence or absence of particular taxa in the phyl- 
losphere was conserved across replicates, suggesting strong selection for community composition, the relative abundance of 
these taxa was highly variable and related to the spatial association of individual plants. Our results suggest that stochastic 
events in early colonization, coupled with dispersal limitation, generated alternate trajectories of bacterial community assembly 
within the context of deterministic selection for community membership. 

IMPORTANCE Commensal bacteria associated with plants help protect their hosts against infection and promote growth. Bacteria 
associated with plant leaves (the "phyllosphere") are highly abundant and diverse communities, but we have very limited infor- 
mation about their ecology. Here, we describe the formation of phyllosphere communities on the plant model organism Arabi- 
dopsis thaliana. We grew a large number of plants in a greenhouse and measured bacterial diversity in the phyllosphere 
throughout the Arabidopsis life cycle. We also measured the diversity of airborne microbes landing on leaves. Our findings show 
that plants develop distinctive phyllosphere bacterial communities drawn from low-abundance air populations, suggesting the 
plant environment is favorable for particular organisms and not others. However, we also found that the relative abundances of 
bacteria in the phyllosphere are determined primarily by the physical proximity of individual plants. This suggests that a mix- 
ture of selective and random forces shapes phyllosphere communities. 



Received 16August2013 Accepted 19 December 2013 Published 21 January 2014 

Citation Maignien L, DeForce EA, Chafee ME, Eren AM, Simmons SL 2014. Ecological succession and stochastic variation in the assembly of /\rafa/dop5/5 thaliana phyllosphere 
communities. mBio 5(1 ):e00682-1 3. doi:1 0.1 1 28/mBio.00682-l 3. 
Editor David Relman, VA Palo Alto Health Care System 

Copyriglit © 2014 Maignien et al.This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-ShareAIIke 3.0 Unported 
license, which permits unrestricted noncommercial use, distribution, and reproduction in any medium, provided the original author and source are credited. 
Address correspondence to Sheri L. Simmons, sherisim@gmail.com. 



Plants in nature are colonized by a large, diverse array of non- 
pathogenic microbes (1). Surveys of both root- and leaf- 
associated microbes demonstrate that these communities are 
highly abundant (2), are species and ecotype specific (2-4), and 
can have significant phenotypic impacts on host plants (5). Here, 
we focus on microbial communities associated with the aerial 
parts of plants, primarily leaves (the "phyllosphere") (2). Phyllo- 
sphere microbes experience high levels of UV exposure, water 
stress, large shifts in temperature (6), and heterogeneous nutrient 
availability (7). Despite this, the global population of phyllosphere 
microbes is estimated to be -10^'' cells (8), and cell densities in the 
phyllosphere are typically around 10^ to 10^ cells/cm^ (6). Phyllo- 
sphere bacteria can protect their hosts against pathogen infection 
(9, 10) and produce plant growth-promoting hormones (11). On 



a larger scale, phyllosphere communities affect biogeochemical 
cycling through the breakdown of plant-released methanol (12) 
and are a major source of bacteria to the atmosphere (13). 

Despite the importance and abundance of phyllosphere mi- 
crobes, little is known about the ecological processes regulating 
their dynamics. Culture-independent surveys of phyllosphere 
bacterial diversity in natural communities have identified seasonal 
variation (14, 15), geographic site (16, 17), and plant species (18) 
as regulators of community structure. In most of these studies, 
however, a significant portion of inter- and intraplant variability 
in phyllosphere composition remains unexplained by environ- 
mental structuring factors (17-19). For example, a study of 
Methylobacterium on natural populations of Medicago and Arabi- 
dopsis found that over half the variance in community structure 
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was not explained by site or plant species (17). These results sug- 
gest that stochastic processes, including dispersal, ecological drift, 
and colonization history, may play a significant role in structuring 
phyllosphere communities. 

The intriguing mix of deterministic (niche-based) and sto- 
chastic processes operating in phyllosphere communities makes 
them interesting model systems to address key questions in mi- 
crobial community assembly (14,20). The balance between deter- 
ministic and stochastic forces in microbial systems is a continuing 
subject of debate (21, 22). Most studies of bacterial assembly in 
natural communities find that both sets of processes are impor- 
tant (23-25) but are limited by low replication and the difficulty of 
controlling assembly in field environments. Laboratory model 
systems, in contrast, allow high levels of replication and experi- 
mental control. They have the potential to generate predictive 
explanations of community assembly (22). 

In this study, we use bacterial communities associated with the 
phyllosphere of Arabidopsis thaliana to explore the repeatability of 
community assembly and the importance of deterministic and 
stochastic processes. The phyllosphere has several advantages as a 
model. Leaf surfaces are complex natural habitats with extensive 
microscale variation (26, 27) and diverse substrates for microbial 
growth (6). We can initialize large numbers of new communities 
on axenic plants and track community assembly through time 
using high-resolution sequencing methods. The emergence of 
new leaves allows us to observe the colonization of new habitats 
and assembly of the resulting microbial communities. Finally, 
phyllosphere communities are easily manipulated and replicated 
in the lab. 

Arabidopsis thaliana, the premier model organism in plant bi- 
ology, is ideal for experimental studies of phyllosphere assembly. 
The availability of inbred Arabidopsis lines reduces the influence 
of host-to-host genetic variability on phyllosphere structure. Ara- 
bidopsis has a short life cycle and can be grown easily in large 
quantities (28). Additionally, it is increasingly used as a model to 
explore plant-microbial community interactions; two recent 
landmark studies oi Arabidopsis rhizosphere microbes discovered 
an influence of host genotype and developmental stage on micro- 
bial community structure (3, 4). 

To test the repeatability of phyllosphere community assembly 
and the influence of deterministic and stochastic forces, we grew 
many replicate plants of the same genotype in a greenhouse envi- 
ronment (Fig. 1). We measured bacterial community structure in 
the phyllosphere and the subset of airborne bacteria passively set- 
tling near leaves throughout the Arabidopsis life cycle. Replicate 
phyllosphere communities converged to a shared community 
composition distinct from air as plants matured, revealing the 
importance of selective forces. However, we also observed that the 
relative abundance of these taxa was clearly related to the spatial 
association of individual plants. Our results suggest a strong in- 
fluence of assembly history and small-scale dispersal limitation on 
phyllosphere community structure. 

RESULTS 

We sampled 32 individual plants and 21 airborne colonizing mi- 
crobial communities for characterization with tag pyrosequenc- 
ing of 16S rRNA genes and used quantitative PGR (qPCR) to 
obtain estimated cell densities in the phyllosphere throughout the 
Arabidopsis life cycle. Environmental conditions were relatively 
constant in the greenhouse throughout the experiment (see 



Fig. SIA in the supplemental material); temperature was held at 
20°C, and humidity fluctuated with outside weather. Day length 
was ~9 h. Plant vegetative biomass continued to increase through- 
out the experiment (see Fig. SIB), with most plants bolting by day 
50. At the end of the experiment, most plants had mature siliques 
but had not yet entered senescence. 

Quantitative PGR. We obtained 16S rRNA operon copy num- 
bers for 26 of 33 plant samples, representing 10 time points (see 
Fig. S2); biomass in the remaining samples was too low for accu- 
rate quantification. Reaction efficiency was 71.6% (_R^ = 0.999), 
possibly due to high primer degeneracy. While the 783F primer 
limited chloroplast amplification, it did not prevent 16S rRNA 
gene amplification from host plant mitochondria. Mitochondrial 
contamination was minimal in sequenced libraries ( < 1 % in 47/60 
datasets, 1 to 5% in 13/60 datasets) and could be discounted in 
copy number calculations. Pyrosequencing libraries with a high 
percentage of chloroplast 16S sequences had lower qPCR copy 
numbers (-R^ = 0.3), suggesting that amplification of chloroplast 
rRNA genes occurred only in the absence of significant bacterial 
abundance. We estimated operon copy number for each opera- 
tional taxonomic unit (OTU) using rrndb (29) and corrected read 
counts for each OTU to generate estimated cell counts in each 
sample (average rrn copy number/sample = 3.6). We observed a 
linear correlation between leaf area and rosette dry weight after 
leaf wash (with trimmed root and stems) across 12 plants (_R^ = 
0.952). We thus used dry weight data to extrapolate leaf area and 
deduce bacterial densities on leaves, taking into account both 
adaxial and abaxial leaf area. 

Microbial population density on leaves was relatively constant 
between days 45 and 67 (1.7 X 10^ to 2.7 X 10* cells per cm=^), 
rising to 4.8 X 10* ± 2.9 X 10" cells/cm^ on day 73 (see Fig. S2). 
Total estimated cell counts per plant rose from 8.0 X 10^ ± 0.54 X 
10^ on day 19 to 2.4 X lO^^ ± 1.3 X 10** on day 73, an ~300-fold 
increase. These results suggest that microbial density on Arabidop- 
sis leaves reached a steady-state density of 10* cells/cm^ roughly 
25 days following initial colonization. This value is at the lower 
end of the reported range for phyllosphere microbial densities 
(30). Because humidity is a strong determinant of cell density in 
the phyllosphere (31), these lower densities are likely due to the 
relatively dry conditions in the greenhouse during our experi- 
ment. 

Verification of sterility. To accurately compare airborne mi- 
crobial immigrants with the phyllosphere community, it was im- 
portant to ensure that contamination was not introduced at any 
point in the sampling and extraction process. To this end, we 
implemented rigorous sterility controls throughout the experi- 
ment. Initially, we verified the sterility of soil and seedlings germi- 
nated from sterile seeds on phytoagar plates. Pyrosequencing of 
the 16S rRNA gene v4v6 region from seedlings extracted prior to 
planting produced >99% organellar sequence (chloroplast and 
mitochondrial 16S rRNA genes). Soil sterility was verified through 
extraction of soil samples prior to planting and amplification with 
universal bacterial 16S rRNA gene primers; visualization of reac- 
tions on a Bioanalyzer high-sensitivity DNA chip confirmed the 
absence of amplifiable DNA. Controls from each sampling and 
extraction event were also amplified with universal bacterial prim- 
ers and visualized on the Bioanalyzer. Controls from each sam- 
pling and extraction event were also amplified with 16S primers, 
and no contamination was observed in any sample. 
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Tag sequencing and diversity analysis. Of the 505,082 raw 454 
GS-FLX reads obtained using primers targeting the v4v6 hyper- 
variable region, 346,866 passed our quality control pipeline and 
were assigned a taxonomic rank using the GAST process (32). We 
excluded sequences not affiliated with bacteria (43,536 reads), in- 
cluding those derived from chloroplasts or mitochondria (0 to 
65%, mostly chloroplast). The remaining sequences, trimmed to 
the v5v6 region, were clustered into OTUs at 97% similarity using 
a usearch6-based clustering algorithm (33). This resulted in an 
abundance matrix of 1,758 OTUs containing 266,043 reads, in- 
cluding singleton and doubleton OTUs. A large number of OTUs 
were unique to either air (1,003 OTUs) or plants (435 OTUs), 
while 320 OTUs were shared by both sample types, suggesting that 
increased sampling depth is required to capture rare OTUs in air 
that colonize plants. A prior comparison of contemporaneous air 
and phyllosphere communities using denaturing gradient gel 
electrophoresis (DGGE), a lower-resolution method, also found 
that few taxa (2/28 bands) were shared between air and leaves ( 19) . 

We calculated within (alpha)-sample diversity on the full ob- 
servation matrix using both parametric (CatchAll) and nonpara- 
metric (Chaol) richness estimators (34). Initial plant community 
richness varied between 79 and 114 at day 19 and increased 6- to 
7-fold by day 60, depending on estimation method (see Fig. S3). 
The estimated richness of the airborne colonizing community 
does not show a time-related evolution and oscillated between 131 
and 314 OTUs. Between (beta) -sample diversity was computed 
using QIIME (35) and Vegan (36) software after abundance nor- 
malization to the minimum sample depth (896 sequences) and 
removal of singleton and doubleton OTUs. The resulting abun- 
dance matrix contained 261,637 sequences and 397 OTUs over 53 
samples, with an average of 4,936 reads per sample. Subsampling 
and removal of singleton/doubleton OTUs did not affect beta di- 
versity results (data not shown). 

Establishment of a "mature" phyllosphere community. We 
analyzed plant and airborne colonizing community structure 
based on OTU presence/absence, using principal coordinate anal- 
ysis (PCoA) on pairwise unweighted beta diversity metrics (Jac- 
card [Fig. 2A] and uwUniFrac [see Fig. S4A]). Airborne coloniz- 
ing communities formed a distinct cluster (Fig. 2 A) separate from 
plant communities and did not show any time-dependent change. 
Phyllosphere communities initially resembled the airborne colo- 
nizing community (day 19) but then departed from the air cluster, 
converging toward a mature phyllosphere cluster containing all 
communities from plants at day 60 and later. The composition of 
day 60+ communities was more phylogeneticaUy homogenous 
than airborne immigrating communities (see Fig. S4A). In addi- 
tion, the laccard distance between replicate plant communities 
decreased with time (Fig. 2B), indicating a convergence of phyllo- 
sphere communities in OTU membership. Permanova analysis on 
Jaccard distance matrices indicated a significant effect of both or- 
igin (air/plant, P < 0.001) and sampling day {P < 0.001) on com- 
munity similarity. Sampling day was highly significant when con- 
sidering plant communities alone {P = 0.005). 

The dominant members of mature phyllosphere communities 
(see Table SIA in the supplemental material) were Acinetobacter 
(29%), Variovorax (12%), Pseudomonas (11%), unidentified 
Sphingobacteriaceae (7%), Rhodococcus, Ochrobactrum, and 
Chryseobacteriutn (4%). The most abundant taxa in air included 
Alicyclobacillus (23%), Ralstonia (26%), Acinetobacter (9%), 
Methylobacterium (5%), and Pseudomonas (3%) (see Fig. S5). We 



identified some important differences between the phyllosphere 
community structures of our greenhouse-grown plants and the 
phyllospheres of field Arabidopsis (2). In particular, we noted an 
underrepresentation of Methylobacterium and Sphingomonas, 
which typically dominate in field populations (37-39), and an 
overrepresentation of Acinetobacter and Pseudomonas taxa in our 
late-stage communities compared to in field plants. 

We identified 51 biomarker OTUs (16 air, 41 plant) that 
strongly differentiated air and mature plant communities using 
LEfSe software (40) (see Fig. S6). Sets of air and plant biomarkers 
are phylogeneticaUy distinctive at the phylum level (see Fig. S6A); 
air biomarkers are primarily in Firrnicutes, while plant biomarkers 
are associated with Bacteriodetes and Proteobacteria. Biomarker 
abundance trajectories confirm that they were consistently more 
abundant in plants or air during the experiment. OTUl (Alicyclo- 
bacillus) dynamics (Fig. 2C) are representative of taxa that were 
abundant in the pool of colonizing microorganisms but that were 
rapidly excluded from the phyllosphere microbiota. Conversely, 
OTU4 (Variovorax) and OTU 10 (Rhodococcus) are representative 
colonists initially rare in air that gradually become dominant 
members of the mature phyllosphere community (see Fig. 2D and 
E). 

We also used LEfSe to identify OTUs appearing at different 
stages of community succession. The number of distinct plant 
biomarkers observed increased as sampling progressed. On each 
sampling day between day 29 and day 50, less than 5 new bio- 
marker OTUs were identified as statistically significant. On day 
60, 19 new biomarker OTUs were identified. These biomarker 
groups represent distinct stages in ecological succession in the 
phyllosphere community (see Table SIB). Notably, some OTUs 
from genera frequently observed in field populations of Arabidop- 
sis (Sphingomonas, Methylobacterium [2]) appeared only later in 
assembly (day 50-1- ). In contrast, several biomarker OTUs in the 
order Sphingobacteriales appeared earlier in community develop- 
ment. 

We also observed several remarkable cases of host discrimina- 
tion between closely related OTUs within the same genus. The 
genus Methylobacterium is one of the most environmentally abun- 
dant genera associated with plants, but the factors regulating its 
distribution are largely unexplained (17). Methylobacterium spe- 
cies are facultative methylotrophs that can grow on a variety of C2, 
C3, and C4 compounds, and some strains produce plant growth- 
promoting hormones (11). We identified strong host selection for 
Methylobacterium OTUs. From a diverse colonizing pool of 13 
OTUs present in air communities, a single one (OTU 31) was 
significantly represented in plants (see Fig. S7). This abundant 
OTU was identified as a biomarker for plant communities (see 
Fig. S7; linear discriminant analysis [IDA] score of 3.4), while 
most others were clearly excluded from Arabidopsis phyllospheres 
in spite of their abundance in air. Indeed, two other Methylobac- 
terium OTUs were identified as air biomarkers, suggesting they 
may have originated from other local plant species and could not 
successfully colonize Arabidopsis. 

Community abundance structure. We examined the struc- 
ture of phyllosphere communities using PCoA on beta diversity 
metrics weighted by OTU abundance (Morisita-Horn [Fig. 3A] 
and wUniFrac [see Fig. S4B in the supplemental material]). The 
structure of airborne colonizing communities, based on Morisita- 
Horn distances, remains relatively stable over time, forming a 
cluster distinct from all but the earliest (day 19) phyllosphere 
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FIG 1 Experimental design. Plants were placed in four growtli trays inside 
two custom-built shade boxes side by side on a greenhouse bench. Arabidopsis 
Col-0 plants were placed in trays 1 and 2 (gray). Trays 3 and 4 (crosses) con- 
tained Arabidopsis gl-1 plants (data not shown). Trays were rotated at every 
watering and sampling event by each tray being moved one place to the right. 
Note that trays 1 and 2 remained in separate shade boxes untU day 55 of the 
experiment. Glass slides coated with CellTak (small black rectangles) were 
placed in the center of each tray and sampled/replaced simultaneously with 
sampling of plants. Flats within each tray were rotated back to front and front 
to back as shown every time trays were rotated. 



communities. This pattern is thus similar to the one observed 
using presence/absence metrics. However, plant communities 
form two groups corresponding to the two trays used for plant 
growth (Fig. 1). Moreover, this tray effect masks the time- 
dependent trajectory apparent in Fig. 2. The distinction between 
trays appears early during community assembly: day 19 plant 
samples are already distinct, and by day 29, plant samples occupy 
a central position within their respective groups (Fig. 3). Per- 
manova analysis identified a significant effect of the tray on 
Morisita-Horn similarity (P < 0.001) but no effect of sampling 
day(P= 0.5). 

LEfSe analysis identified 22 OTUs as biomarkers for either tray 



1 or 2, with a prominent role of Pseudomonas in tray 1 and Acin- 
etobacter in tray 2 (Fig. 3). OTUs from both genera were among 
the strongest tray biomarkers identified. These taxa strongly dom- 
inated phyllosphere communities; individual OTUs from these 
genera appeared at relative frequencies of up to 52% (Pseudomo- 
nas, day 50) and 72% (Acinetobacter, day 24). Notably, Pseudomo- 
nas and Acinetobacter OTUs were very abundant early in commu- 
nity development (Fig. 3B and E). LEfSe also identified known 
plant- and soil-associated taxa as tray-specific biomarkers, includ- 
ing Rhodococcus, Methylobacterium, Novosphingobium, and 
Chryseobacterium (Fig. 3C and D; see also Fig. S8 in the supple- 
mental material), which appear later in the time series. 

Ecological differentiation of phyllosphere bacteria at the 
sub-OTU level. Interestingly, the strong tray pattern we observed 
using the Morisita-Horn diversity index almost completely col- 
lapsed when the phylogeny-based wUniFrac distance was used 
(see Fig. S4B). This shows that bacteria accounting for differences 
between trays are closely related phylotypes. For example, the 
dominant OTUs driving community abundance patterns belong 
to the genera. Acinetobacter and Pseudomonas, which are both part 
of the same order (Pseudomonadales) . 

We hypothesized that these closely related taxa may be ecolog- 
ically equivalent (41) in phyllosphere communities. We examined 
whether ecological processes structuring bacterial communities 
also operate at the sub-OTU level (corresponding to related spe- 
cies, strains, or genotypes of the same species). OTU clustering, 
while minimizing the effect of sequencing errors on microbial 
diversity analysis, can obscure biologically significant variations 
between organisms with highly similar (>97%) 16S rRNA gene 
sequences. We thus used oligotyping, a powerful new method, to 
better distinguish true biological variations from sequencing error 
noise (42). Oligotyping uses Shannon entropy analysis to isolate 
meaningful positional variation within sets of similar 16S rRNA 
gene sequences that with standard methods would converge into 
the same cluster; each unique set of variants is identified as an 
ohgotype. 

The 20 most abundant genera, which together accounted for 
99 OTUs and 77% of reads in the data set, were individually de- 
composed by this method, identifying 234 different oligotypes. 
For six of these genera, oligotype decomposition matched OTU 
clustering. In all other cases, however, oligotypes demonstrated 
distinct distributional patterns. As an example, Acinetobacter 
OTU 0, the most abundant plant-associated taxon overall, could 
be decomposed into four oligotypes that each exhibited a dramat- 
ically different spatial distribution (Fig. 4) . In spite of the presence 
of oligotypes 1 and 2 in airborne immigrating communities and in 
initial phyllosphere communities, only oligotype 1 successfully 
colonized the tray 2 phyllosphere. These results show that two 
Adnetofoflcfer phylotypes (identified as Acinetobacter johnsonii and 
Acinetobacter junii based on BLAST results) with 98.2% sequence 
similarity in the 16S rRNA gene had opposite colonization pat- 
terns. 

DISCUSSION 

In this study, we characterized the development of bacterial com- 
munities in the phyllosphere of Arabidopsis thaliana over the plant 
life cycle. We used a large number of replicate plants grown in a 
controlled environment to test the repeatability of community 
assembly, the first such study of phyllosphere assembly to date. 
Our results demonstrate that the presence or absence of particular 
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FIG 2 Mature plant-associated microbial communities have a distinct membership from air communities and show a clear trajectory over time. (A) Principal 
coordinate plots using the membership-based Jaccard index to measure beta diversity, based on OTUs clustered at 97% similarity. Each dot represents a single 
community (blue, air; green, plant). Dot size is scaled by sampling day (day 19, small; day 73, large). Arrow indicates the trajectory from early to mature 
phyllosphere communities. (B) Mean Jaccard distance between replicate plants on each sampling day shows increasing similarity with time. (C) Abundance 
trajectories over time of the day 60-1- air biomarker OTU 1. Blue indicates relative abundance (percentage of total read count) in air, and green indicates relative 
abundance in plants. Error bars show the standard deviation across triplicates (plants) and duplicates (air). Abundance trajectories over time of day 60-1- plant 
biomarker OTU 4, a day 55 biomarker (D), and OTU 10, a day 50 biomarker (E). Coloring is identical to that described for panel C. 



taxa in the phyllosphere is under strong deterministic selection 
(Fig. 2). The membership of phyllosphere communities initially 
mirrored immigrant airborne microbes but subsequently con- 
verged to a phylogenetically distinctive community composition. 
We observed strong, reproducible successional dynamics in com- 
munity membership. In contrast, the relative abundance of bac- 
terial taxa in the phyllosphere was highly variable among repli- 
cates and was strongly related to the spatial association of 
individual plants (Fig. 3). Separation by spatial association began 
very early in community formation and continued throughout the 
plant life cycle, generating alternate trajectories of community de- 
velopment. These results suggest that stochastic forces play a sub- 
stantial role in structuring phyllosphere communities. 

Stochastic forces in community assembly. A recent concep- 
tual synthesis (43) identified four main processes in community 
assembly: selection, drift, speciation, and dispersal. Only one of 
these forces, selection, predicts a correlation between habitat 
niche structure and community composition. The other three sto- 
chastic processes can operate in conjunction with selection or to- 
gether to generate the classical neutral model of Hubbell (44). 
Theory suggests that stochastic variation in colonization order can 
have a significant impact on community assembly, resulting in 
high beta diversity among similar sites (45, 46). Dispersal limita- 



tion can reinforce the effects of colonization history on beta diver- 
sity (46). Implicit in these models is the assumption that niche- 
based selection operates on sets of functionally equivalent 
colonizers to generate alternate community structures in similar 
sites. Stochastic niche theory (47), for example, predicts that pri- 
mary colonizers rapidly occupy the broadest available niches and 
preempt further colonization by similar species. Secondary colo- 
nizers must be able to fit into the remaining niche space and grow 
rapidly enough to overcome drift. Under this model, alternate 
states are easily generated given a sufficient diversity of primary 
colonizers. Alternatively, stochastic variation in colonization or- 
der could lead to priority effects, where the chance arrival of a 
particular early colonizer alters the habitat in a way that favors the 
growth of specific secondary colonizers (48). 

Our observations suggest that stochastic colonization dynam- 
ics and dispersal limitation played a central role in shaping the 
abundance structure of phyllosphere bacterial populations. The 
convergence in community membership across replicate plants 
over time indicates that host-microbe and/or microbe-microbe 
interactions combined to shape niches favoring the growth of par- 
ticular taxa. It further suggests that replicate plants have similar 
niche structures. Given this result, it is difficult to interpret the 
strong effect of spatial association on bacterial abundance (Fig. 3) 
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FIG 3 Spatially associated plants share similar taxon abundances. (A) Principal coordinates plot using the abundance-based Morisita-Horn index to measure 
beta diversity, based on OTUs clustered at 97% similarity. Blue, air; green, tray 1 plants; orange, tray 2 plants. Scaling is the same as described for Fig. 1. Arrows 
indicate trajectories from newly colonized plants to tray-specific communities. (B) Abundance trajectory of the dominant tray 1 biomarker Pseudomonas; error 
bars indicate the standard deviation of replicate plants. Green bars indicate relative abundance of the OTU in tray 1, orange bars in tray 2, and blue bars in air. 
(C, D) Abundance profiles for major tray marker taxa identified by LEfSe at successive time points. (C) Abundance of the tray 1 biomarker Rhodococcus. (D) 
Abundance of the tray 2 biomarker Methylobacterium. (E) Abundance of the dominant tray 2 biomarker Acinetobacter. Note that three plants were sampled at 
every time point, chosen randomly; possible configurations for sampling between trays 1 and 2 were (2,1), (1,2), (3,0), or (0,3). 



as resulting from purely niche-based forces. The two experimental 
trays were treated identically and frequently rotated, yet plants in 
one tray had phyllosphere communities dominated by Acineto- 
bacter while communities in the other tray were dominated by 
Pseudomonas. These two highly abundant taxa appeared very early 
in community assembly (Fig. 3B and E). They are the strongest, 



and earliest appearing, tray biomarkers identified by LEfSe (see 
Fig. S8). Distinct tray community structures were already appar- 
ent at 5 days postcolonization (Fig. 3A) and continued until the 
end of the experiment. 

The most parsimonious explanation is that random initial col- 
onization events led to the early dominance of each taxon in a 
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particular tray, which is supported by the divergence of tray abun- 
dance structures immediately after colonization (day 19; Fig. 3). 
Populations developing from these early events were maintained 
on the initial leaves and then propagated to newly formed leaves. 
Given that trays were separated by roughly 1 m, and plants within 
a tray were separated by at most 40 cm, it is plausible that dispersal 
limitation reinforced alternate community trajectories through- 
out the experiment. This hypothesis is supported by the tray- 
specific distribution of several late-appearing colonizers, includ- 
ing OTUs from known plant- and soil-associated genera 
{Rhodococcus, Methylobacterium, Novosphingobium, and Chryseo- 
bacterium; Fig. 3C and D). These observations also fit an impor- 
tant prediction of the stochastic niche model (47): initial coloniz- 
ers have access to more resources and are more likely to develop 
dominant populations compared to competitors with smaller 
niche utilization capabOities. 

The alternative assembly trajectories in trays 1 and 2 may rep- 
resent an interesting case of ecological drift between spatially sep- 
arated communities (22). It shows remarkable similarity to other 
cases of stochastic community assembly, such as the dependence 
of mouse gut microbiome composition on housing cage (49) . Dis- 
persal limitation as a key factor in phyllosphere assembly is also 
consistent with field surveys identifying a large effect of geo- 
graphic site on phyllosphere structure (17). An alternate possibil- 
ity, which does not require the assumption of dispersal limitation, 
is that Pseudomonas and Acinetobacter each altered the leaf surface 
habitat in such a way to favor the growth of particular sets of 
secondary taxa. 

Functional equivalence in microbial populations. Our obser- 
vation of alternate phyllosphere assembly trajectories in spatially 
separated plants raises the obvious question of whether these 
states are functionally equivalent. Empirical evidence from other 
microbial systems is mixed. A combination of inoculum source 
and housing cage led to differences in mouse gut microbiome 
function (49), and stochastic colonization coupled with biotic in- 
teractions led to differences in function across artificial microbial 
reactors (50). On the other hand, there are many examples of 
variation in bacterial community structure across similar sites 
(e.g., multiple wastewater reactors) coupled with conservation of 
function (23, 51). 

The functional equivalence of cooccurring taxa is central to the 
neutral theory of community assembly (52). Hubbell (53) argued 
that sets of functionally equivalent species can evolve within the 
niches that are most prevalent over evolutionary time. Therefore, 
the dynamics within each niche are neutral, but selection has oc- 
curred to set the boundaries of each species set. The necessary 
condition for this to occur is the absence of factors that promote 
competitive exclusion between functionally similar species. Other 
theoretical work suggests that this pattern of sets of functionally 
similar species ("emergent neutrality") can appear via several dif- 
ferent possible pathways (54, 55) and is more likely to occur in 
species-rich communities (56). 

While our experiment did not directly test for functional 
equivalence, some of our findings are suggestive. The dominant 
early colonizing tray biomarker taxa, Acinetobacter and Pseu- 
domonas, arguably fill similar niche spaces in the leaf environment 
as primary colonizers. Acinetobacter and Pseudomonas have on 
average 6 and 5 ribosomal operons per genome (29), respectively, 
making them classic r-strategist taxa capable of rapid growth in 
response to nutrient availability (57). Additionally, species within 



both genera are frequently found in association with plants and 
have growth-promoting properties (58, 59). However, our analy- 
ses also identified examples of ecological selection between closely 
related taxa. Two oligotypes of Acinetobacter, which have 98.2% 
similar 1 6S rRNA gene sequences, exhibited distinct distributional 
patterns; while both were present in the airborne immigrating 
community, only one was able to establish on tray 2 plants (Fig. 4) . 
We also found that when multiple OTUs of a particular genus 
were present in airborne immigrants, often only one successfully 
established on plants (e.g., Methylobacterium; see Fig. S7 in the 
supplemental material). Testing the hypothesis of functional 
equivalence between alternate states will require additional repli- 
cated experiments, as well as direct assays of community function 
(for example, through shotgun metagenomics) and community- 
level competition experiments. 

Conclusion. Through detailed characterization of 
Arabidopsis-associated leaf microbiota and airborne colonizing 
microbes over the plant life cycle, we identified key ecological 
forces driving microbiome assembly. On the one hand, conver- 
gence in microbiome membership as plants mature indicates that 
plants exert a strong selective force on the identity of colonizing 
microbes (who can colonize). However, variation in the abun- 
dance — as opposed to the presence — of dominant taxa is strongly 
related to spatial associations between plants. This variance is best 
explained by stochasticity in initial colonization events and sub- 
sequent limited dispersal, as predicted by different community 
assembly models (neutral or stochastic niche). Further experi- 
mentation with controlled community assembly is required to 
assess the repeatability and robustness of these coupled niche and 
stochastic dynamics. 

Our results demonstrating an interaction of niche and stochas- 
tic effects are suggestive, but it will be necessary to replicate our 
experiment in different controlled environments (research green- 
houses or hoophouses) and with a larger number of plants to test 
their generality with respect to microbial colonization of the phyl- 
losphere. Our findings of two alternate community structures also 
raise the question of how many community trajectories are possi- 
ble in the phyllosphere environment and whether they differ in 
relative fitness. Again, similar experiments with a larger number of 
replicate plants will be needed to address these questions. Addi- 
tionally, the conspicuous role of dispersal patterns suggested by 
our results needs to be directly tested and quantified. We antici- 
pate that more complex experimental designs using various dis- 
persal rates between phyllosphere metacommunities will bring 
valuable insights to this problem. 

Our results have implications for the design of experiments 
aiming at testing the effects of particular treatments or environ- 
mental variables on phyllosphere community structure. We show 
that both stochastic events and dispersal limitation can account 
for significant beta diversity between spatially separated replicate 
communities. Large pools of replicates are necessary to account 
for this inherent stochasticity, and randomization of control and 
tested individual plants is necessary to avoid confounding ecolog- 
ical drift among pooled individuals with experimental treatments. 

Overall, this study provides a novel ecological framework for 
studies of microbiome assembly and points to the importance of 
highly replicated, controlled longitudinal studies of microbial 
community development. 



January/February 2014 Volumes Issue 1 e00682-13 



mfiio' mbio.asm.org 7 



Maignien et al. 



MATERIALS AND METHODS 

Experimental design. We initiated a 73-day time series experiment using 
72 microbe-free Arabidopsis thaliana Col-0 plants planted in sterile soil in 
a greenhouse (Fig. 1). Every 4 to 6 days, three plants were destructively 
sampled for phyllosphere community analysis. In addition, we used glass 
slides coated with adhesive protein (CellTak) in order to sample the air- 
borne community. The slides were located among the plants at the level of 
plant leaves. The slides were intended to act as passive traps for microbes 
arriving on the leaf surface during a particular time interval rather than be 
representative of the total airborne microbiota. They allow us to deter- 
mine, to some extent, the composition of the immigrant community in 
the absence of the dynamic microbe-microbe and microbe-host interac- 
tions shaping leaf surface communities. 

Plant germination and growth. Seeds of Arabidopsis thaliana were 
surface sterilized, germinated on standard phytoagar plates, and trans- 
ferred after 14 days of growth to sterile soil in a Conviron greenhouse 
(Falmouth Technology Park, Falmouth, MA). The light regime was 
roughly 9 h of light and 1 5 h of dark. Sterile soil was obtained by saturating 
dry soil (Lehle Seeds, Round Rocks, TX) with water in an autoclave bag. 
After 48 h at room temperature, the wet soil was autoclaved twice for 
45 min at 24-h intervals (60). Autoclaved soil was tested for the presence 
of amplifiable bacterial 16S rRNA genes. Prior to being planted, seedlings 
were sampled from phytoagar plates and tested for sterility with 1 6S rRNA 
gene PGR. Upon planting in soil, each plant was transferred to a separate 
well of a 6-well insert. Inserts were placed in plastic trays under shade 
boxes constructed with wire screens to reduce light levels to those stan- 
dard iox Arabidopsis growth (28). The first sampling occurred on day 19, 
5 days after transfer from plates to soil. 

Environmental conditions in the shade boxes were monitored using 
HOBO sensors (temperature/RH [relative humidity] and light as PAR 
[photosyntheticaUy active radiation]) and a MicroStation data logger 
(Onset Computer, Bourne, MA). The greenhouse temperature set point 
was 20°C, and humidity was not controlled. Plants were watered every 
3 days and fertilized with 200 ppm 20/10/20 fertilizer IX per week starting 
at day 27. Plant and tray positions were randomized at each watering, 
fertilization, and sampling event. Plants were randomized within each 
tray, and trays were moved between the two shade boxes, but plants were 
not moved from tray to tray (Fig. 1). 

Air and phyllosphere sampling. We captured airborne microbes rep- 
resenting the colonizing community using sterilized standard microscope 
slides coated with the biological adhesive CellTak (BD Biosciences) at a 
density of 3.5 p,g/cm^. We constructed small platforms to hold slides at the 
level of plant leaves. Slides were left in place during each sampling interval 
(~5 days), collected, and replaced at the next sampling event. In prelimi- 
nary experiments, fluorescent microscopy using 4',6-diamidino-2- 
phenylindole (DAPI) staining indicated the presence of intact cells on 
slides exposed to air for ~l-week periods (data not shown). Slides were 
incubated vnth trypsin to remove CellTak and detach microbes; the wash 
solution was collected on sterile 0.22-fj,m-pore-size filters, and DNA was 
extracted using the Biostic bacteremia DNA extraction kit (MoBio Labs). 
Due to the low biomass present in air, we used rigorous sterile technique 
and extensive negative controls to exclude contaminant microbes (see 
Results). 

Three whole Col-0 plants, randomly selected, were sampled at each of 
1 1 time points. Plants were randomly drawn from both trays, resulting in 
an unbalanced number of plants from each tray at different time points. 
Roots and flowering stems were trimmed off, and then rosettes were 
placed into a solution of 0.2% Silwet in TE (10 mM Tris, 1 mM EDTA, 
pH 7) (37, 61 ) and sonicated in a bath sonicator for 10 min. Rosettes were 
removed, dried at 70°C overnight, and subsequently weighed. Wash solu- 
tions were preflltered through 5-/L(,m-pore-size filters and then collected 
on 0.22-p,m-pore-size sterile filters. DNA was extracted from filters with 
Biostic bacteremia kits. Plant leaf area was calculated based on dissected 
plant photos using Image) software (62). Plant dry weight was measured 
on a high-precision scale after overnight drying of the plants at 70°C. It 



was not possible to obtain wet weights due to the addition of Silwet/buffer 
to plants in the wash procedure. One of the day 60 plants and aU day 45 
CellTak slide replicates were omitted after manipulation problems 
(breakage of the filtration membrane). 

Sterility tests. During each sampling event, we collected four negative 
controls to detect possible contamination during plant and air sampling. 
First, we filtered sterile Silwet wash solution through our Swinnex appa- 
ratus onto a 0.22-/j,m-pore-size filter. Second, we sampled 500 of Silwet 
wash solution alone. Third, we filtered 30 ml of the trypsin solution used 
in CellTak slide extraction onto a 0.22-f.im-pore-size filter. Finally, prior 
to each sampling event, we placed a sterile CellTak slide inside a 50-ml 
Eppendorf tube within the laminar flow hood. This control slide was 
transported to the greenhouse, placed on the bench near exposed slides, 
and subsequently analyzed. We also included a blank control filter in every 
batch of DNA extractions, as well as standard PCR-negative controls. All 
filters were extracted with MoBio bacteremia DNA kits, and the resulting 
DNA was amplified with standard v4v6 primers. The lack of amplification 
signals from DNA extracted from sterilized soil using a Bioanalyzer high- 
sensitivity DNA chip (detection limit of 5 pg//J.I, manufacturer's specs) 
was interpreted as sufficient evidence for soil sterility. Because sterilized 
soil exposed to air for 5 days showed a strong PCR product using the same 
assay, we conclude that our initial results are not due to inhibition of PCR 
amplification by elements from soil coextracted with DNA. 

qPCR. A SYBR green quantitative PCR (qPCR) assay was developed 
for quantification of bacterial 16S rRNA gene copy number among rosette 
leaf washes using a modifiedl6S rRNA gene primer that limits chloroplast 
16S amplification coupled with the universal primer 1046R (63). A 10- 
fold serially diluted, 5-point standard curve (range, 3 X 10' to 3 X 10^) 
was generated with plasmids containing 16S rRNA gene 783F/1046R in- 
serts from 9 microbial species of various GC content. The standard curve, 
environmental samples, and NTCs (no template controls) were run in 
triplicate on a StepOnePlus real-time PCR system (Life Technologies). All 
rosette samples were run on the same qPCR plate to avoid variation be- 
tween assays due to variable efficiency. Representative sequences for each 
operational taxonomic unit were submitted to rrndb (29) to obtain a copy 
number correction factor to translate qPCR counts into estimated cell 
number. 

Amplification and 454 sequencing of the v4v6 region of 16S rRNA 
genes. The v4v6 variable region of 16S rRNA genes was amplified in trip- 
licate with 454 fusion primers containing adapters and bar codes with 
bacterial primer sequences 518F (5' CCAGCAGCYGCGGTAAN 3') and 
1046R (5' GGAGRRCCATGCANCACCT 3') and sequenced on a 454 
GS-Ti instrument at the MBL as previously described (64). Individual 
plant samples were separately amplified, but CellTak slides from each 
shade box were pooled prior to amplification due to low biomass. 

Bioinformatic analyses. Processing and filtering of v4v6 pyrose- 
quencing reads were carried out using a standard MBL pipeline (32, 64). 
Reads were trimmed to the v5v6 region using a conserved anchor se- 
quence due to low quality at the v4 end, and aU subsequent analyses were 
performed on these data set. Sequences were clustered using a standard 
Usearch6 pipeline into 97% similarity OTUs (33). We used the Catchall 
software (34) for both parametric and nonparametric richness estimation. 
We then used QIIME (35) and Vegan (36) software for beta diversity 
analysis based on a corrected abundance matrix: (i) OTUs containing 1 or 
2 reads were discarded in order to further reduce the number of spurious 
OTUs generated by errors introduced during PCR and sequencing, and 
(ii) observation counts were subsampled to the number of reads present in 
the smallest library (895 reads) for calculation of beta diversity indices. We 
used LEfSE software (40) to identify biomarker OTUs. We used oligotyp- 
ing (42) to identify significant sub-OTU-level variation in each of the 20 
most abundant genera and analyzed the abundance pattern of oligotypes 
across sample sets. 

Nucleotide sequence accession number. AU 16S rRNA gene se- 
quences described in this study have been deposited in the VAMPS ar- 
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chive (http://vamps.mbl.edu) with the accession number 
"SLS_PHY_Bv6v4." 
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